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Abstract. With the aim of understanding the cosmic velocity fields at large scale, we 
investigate the dynamics of a pressureless distribution of gravitational sources moving 
under an anisotropic generalization of Hubble expansion and constraint by Euler- 
Poisson equations system. As a result, it turns out that such a behavior requires the 
distribution to be homogenous, similarly to Hubble law. Among several solutions, we 
show a planar kinematics with constant (eternal) and rotational distortion, where the 
velocity field is not potential. Within this class, the one with no rotational distortion 
identifies to a bulk flow. To apply this model within cosmic structures as the Local 
Super Cluster, the solutions are interpreted as approximations providing us with an 
hint on the behavior of the cosmic flow just after decoupling era up to present date. 
Such a result suggests that the observed bulk flow may not be due exclusively to tidal 
forces but has a primordial origin. 
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1. Introduction 



In the past, the investigations of cosmic velocity fields from redshifts surveys were not as 
successful as that for providing us with powerful tools for the estimation of cosmological 
parameters and the local density fields, see e.g. p. Later on, least action methods 
|2j were also used in reconstruction procedures with more reliable estimates, see e.g. 
[HI IH El U\- Both approaches assume that the peculiar velocity field decomposes 
into a "divergent" component due to density fluctuations inside the surveyed volume, 
and a tidal (shear) component, consisting of a bulk velocity and higher moments, due 
to the matter distribution outside the surveyed volume [Hj. These hypotheses are 
justified by the properties that the irrotational linear perturbations dominate with time 
EH] together with Kelvin theorem which ensures the irrotational characteristics of 
motions. However, it turns out that theses investigations might provide us with biased 
results because of their inherent dependence of sampling characteristics. For example, 
the presence of an high region (Great Attractor) was proposed for accounting of an 
unexpected feature in the cosmic velocities [TIJ^JEIIE], wnen from recent surveys 
Shapley Concentration seems to be a better candidate, although it is not responsible 
for all of the SMAC flow fSl EH E|- Since such an interpretation is not as clear 
cut as that, one might ask whether the presence of a bulk flow originates necessarily 
from a density excess in the spatial distribution of gravitational structures or another 
alternative could be envisaged. For answering this question, we investigate Euler Poisson 
equations system describing an anisotropic Hubble flow of pressureless distribution of 
gravitational sources in a Newtonian schema. 



2. Dynamics of the cosmic expansion 

The dynamics of a pressureless distribution of gravitational sources (dust) is investigated 
by assuming that the motion of sources satisfy the following kinematics 

f=Af , A(t ) = l (1) 

where A = A(t) stands for a 3x3 matrix depending on cosmic time t which has to be 
determined by an observer at rest with respect to Cosmological Background Radiation 
(CMB). The present investigation limits on collisionless motions, which is ensured by a 
non vanishing determinant of matrix A, and because it reads as a unit matrix at given 
t = t a one has det A > 0. Hence, the velocity field v = v(r,t) is given by 

v= — = AA~ 1 f, A = —A (2) 
dt dt w 

where A -1 stands for the inverse matrix (A -1 A = AA _1 = 11). Moreover, accordingly 

to Hubble law, we assume a radial acceleration field 



H.H. Fliche et al. 3 

We assume now that these motions are constraint by Euler equations system 

%= -div (fX) (4) 
dv dv _ _ . , 

where p = p(f, t) stands for the density. By using the trace, the determinant and the 
equality ^detA = Tr f AA _1 j detA, eq. (JH) transforms 



dp -> — 
p = — + v grad p 
at 



= - p div (AA-V) = -p tr (AA- 1 ) = -3p - (6) 
where the doted variables stand for time derivatives, and 

a(t) = vMetA, a Q = a(t ) = 1 (7) 
for the (generalized) expansion factor. Hence, with eq. (JJ) one obtains 

^ r ' *) = § = T > ^° = P° r ° 8 

which simply accounts for the mass conservation p dr^ = pdr 3 = pdetAdr^. 
One has a unique matrix decomposition 

H h 

AA _1 = HI H ^B, #=~, H = H(t ), tr B = (9) 

a z a 

where H = H(t) acts as the usual Hubble factor and B = B(£) stands for a traceless 
matrix herein called distortion matrix. It characterizes a deviation from isotropy of the 
(dimensionless) velocity field, its amplitude is defined by the matrix norm 



||B|| = ytr(B*B) (10) 
where the sign stands for the matrix transposition. For convenience, let us write 

P n (t) =tr(B n ), n = l,2,3 (11) 
According to eq. ()2I5I9|) . one has 

? =((M) 1+ 5(B + §^)) ? (12) 

We assume that the gravitational field satisfies Poisson-Newton equations 

divg = - AixGp + A (13) 

rot# = (14) 

where G is Newton constant of gravitation and A the cosmological constant. According 
to eq. (JI2J), the left hand term of eq. (fT3Jl reads 

di ^=^(f)=3(^ + ^ 2 ) + ^ (15) 

which does not depend on spatial coordinates. Hence, eq. (|13|) tells us that the space 
distribution of sources is homogenous, i.e. 

P = Pit) (16) 
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where 

fio = fS Po ' Ao = 3^' H o = H (to) (18) 

are motion parameters. By multiplying each term of eq. ()17j) by 2aa one easily identifies 
the following constant of motion 



a ° i?? 3 A a 3 



K = h A a - — - - / —da = fi Q + A D - 1 (19) 



Hence, the chronology is given by 



1 oda . . 



(22) 



where H Q > accounts for an expansion (according to observations, the case H < 
which accounts for a collapse is not envisaged), and 

P(a) = A D a 4 - n a 2 + tt Q a - -a 2 f* ^|da > 0, P(l) = 1 (21) 

3 Ji cr 

The constraint given by eq. ()14j) can be written in matrix form as follows 

dg__ fdgV 

dr \df I 

Such a symmetric property with eq. (|12I17|) shows that the matrix 

B = B 2 -^ 2 11+^B (23) 

is traceless and symmetric 

trB = 0, B* = B (24) 
According to eq. 1)3112123124}) . since the field g, which reads 

is radial, the matrix B must be scalar and because it is traceless, one has necessarily 
B = (26) 

2.1. Reference map 

Instead of (t,r), it is more convenient to analyze the dynamics of the cosmic flow in the 
(r, q) coordinates defined by 

dt 

dr = H - (27) 
a z 
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herein called reference map. According to eq. (|2I2(JI23I26)) . the equations of motion read 

da , . 

dr = — ; (30) 



where the distortion matrix B satisfies 

The resolution of these equations can be performed by mean of numerical techniques; 
having solved eq. (j31|) . which gives the evolution with time of distortion matrix B, the 
particles trajectories r i— > <f(r) are obtained by integrating eq. (J29|) and the evolution of 
the generalized expansion factor a from eq. (|3U|). 



2.2. Discussion 

As a result, it is interesting to mention that (as it is the case for Hubble law) this 
anisotropic generalization accounts for homogeneous space distributions of matter. 
Hence, one is forced to ask whether it describes correctly the dynamics of cosmic 
structures because of the presence of strong density inhomogeneities in the space 
distribution of galaxies catalogs. In principle, such a remark should be also sensible 
to question Hubble law when, regardless the isotropy, it is a fact that perturbations are 
not so dominant otherwise it would never have been highlighted. Actually, homogeneity 
is implicitly assumed for the interpretation of CMB isotropy and the redshift of distant 
sources, which provides us with an expanding background. Namely the comoving 
space of FL world model onto which the gravitational instability theory is applied for 
understanding the formation of cosmic structures. It is with such a schema in mind that 
this anisotropic Hubble law may provides us with an hint on the behavior of the cosmic 
flow from decoupling era up to present date in order to answer whether the observed 
bulk flow is due exclusively to tidal forces. 



3. Analysis of analytic solutions 

In this section, we investigate some analytic solutions of eq. (j29l3QI31|) that are obtained 
thanks to particular properties of distortion matrix B. The parameters A , fi and 
H Q given in eq. (jl8j) correspond to cosmological Friedmann-Lemaitre (FL) world model 
parameters. Moreover, the constraint (3% = in eq. (j2()J21j) provides us with the FL 
chronology, where k q given in eq. (JT9j) represents the curvature parameter in the FL 
model (i.e. the dimensionless scalar curvature Q%. of the comoving space, see [IB])) while 
the flatness of (simultaneous events) Newton space. It must be noted that the particle 
position q as defined in eq. ()28|) does not identify to the usual FL comoving coordinate 
because the (generalized) expansion factor a depends on the anisotropy unless fa — 0. 
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3.1. Evolution of functions P n= 2,3 

Because B is a traceless matrix, its characteristic polynomial reads 

Q{s) = det {st — B) = s 3 — l -f3 2 s - ^ 3 (32) 

according to Leverrier-Souriau's algorithm (2D]- With Cayley-Hamilton's theorem {i.e. 
Q(B) = 0) and eq. (}3*T)) we obtain the following differential equations system 

-^2 = - 2(3 3 (33) 
dr 

lh = " 5* (34) 
and we note that the discriminant of third order polynomial Q, it is proportional to 

a = 3$ - ±P! (35) 

is a constant of motion (i.e., da/dr = 0). The integration of eq. (|33I34|) is performed by 
defining p 2 by a quadrature 

r = t + e / . =, e = ±1 36 

and hence (3% from eq. (jUID); in addition of the singular solution 

/?2 = & = 0, {i.e., B 3 = 0) (37) 
defined equivalently either by @ 2 — or /5 3 = 0, according to eq. 



The related dynamics depends on roots ^=1,2,3 of characteristic polynomial Q given 
in eq. (|32jl . ie. the eigenvalues of distortion matrix B. Their real values identify to 
dilatation rates at time r toward the corresponding (time dependent) eigenvectors (not 
necessarily orthogonal). Their sum is null {j3\ = 0) and their product (/% = 3detB) 
is either decreasing with time or is null, according to eq. (|34|). The sign of a given in 
eq. (J35|) is used to classify the solutions as follows : 

• if a = then Q has a real double root r\\ = 772 and a simple one 773. The related 
instantaneous kinematic shows a planar-axial symmetry (either a contraction within 
a plane with an expansion toward a transverse direction or vice versa), see sec. 13.31 
If 771 = 773 then both vanish and the related solution identifies to the singular one 
defined in eq. Q37|) . see sec 13. 2 1 

• if a > then Q has a single real root rji, 

• if a < then Q has three distinct real roots 7/1=1.2,3 • Their order is conserved during 
the evolution (since a coincidence of eigenvalues makes a = 0), the largest one must 
be positive while the smallest one must be negative (because Pi =0). 
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3.2. Planar kinematics 

The singular solution B 3 = shows a FL chronology and the distortion matrix 

B = -B 2 r + B , B 3 = (38) 

is solely defined by its initial value B , according to eq. (|31|). It is neither symmetric 
nor asymmetric (otherwise it vanishes), see eq. Hence, eq. ()29j) transforms 

f r = (-B 2 D r + Bo) q (39) 
which accounts for eternal motions 

q = exp ("B^y + B r^ = (1 + B r) <f (40) 

The trajectory of a particle located at initial position q Q identifies to a straight line 
toward the direction B q . The analysis of B D range (i.e., its image) provides us with 
characteristics of trajectories flow. The nilpotent property of B D shows that its kernel is 
not empty Ker(B Q ) ^ 0. Its dimension dim (Ker(B Q )) = m characterizes the kinematics, 
which is either planar (m = 1) or directional (m = 2), i.e. a bulk flow. Conversely, if 
the kernel of distorsion matrix B is not empty then /5 3 = 3 det (B) = 0, and thus /?2 = 0, 
see to eq. ()33p34|) . Therefore, all planar kinematics can be described by such a model. 

3.3. Planar- Axial kinematics 

If fa 7^ then the chronology differentiates from FL one. Let us focus on the a = with 
two distinct eigenvalues rji ^ 773 class of solutions. With eq. (|HEJ), eq. ()33p34j) integrate 

Pn = -, rn, (n = 2,3), n = r + v/6/3 2 - 1/2 (t ) (41) 

\ T - r J 

which shows a singularity at date r = > that splits the motion in two regimes 
r < r* and r > r*. The complete investigation of this singularity problem demands to 
solve an integro-differential equation, see eq. (j2QI21|) . The roots of Q read 

(r* - r) 

where 771 stands for the double root. Among others, two class of solutions are defined 
by mean of a constant (time independent) matrix P, the projector associated to rji, 
see [201, 

P 2 = P, trP = 2 (43) 
They describe distinct kinematics depending on whether matrix B is diagonalizable. 
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3.3.1. (Irrotational) motions If B is diagonalizable then 

B = 771 (3P - 211) (44) 
From eq. ()29|44|) . one has 

p- = Vl (3P-21)q (45) 
dr 

and the solution reads 

q=- m P^+\(l-P)^ (46) 

m 

where £ is constant. If the eigenvectors are orthogonal then the kinematics accounts for 
irrotational motions. 

3.3.2. Rotational motions If B is not diagonalizable then 

B = rjt (3P — 21) + -rN (47) 

m 



Vl (3P — 211) + — N q (4? 



dr V m 
where N is a constant nilpotent matrix, which accounts for rotational motions on the 
eigenplane of P. The solution reads 

q = -77iPf* + 4 1 - P ) £ + 4 N £ ( 49 ) 

where £ is constant. 

4. Application to flat large scale structures (B 3 = 0) 

The B 3 = class of solutions has interesting properties with regard to the stability of 
large scale structures that show a flat spatial distribution. To answer the question of 
whether observations define unambiguously the kinematics, the distortion matrix B is 
decomposed as follows 

B = S+j(w), trS = (50) 

where S and j(u;) stand for its symmetric and its asymmetric^ component, and 

a 2 ~ -. 1— - 
oj = — a, a = - rot v (51) 

H 2 



accounts for the motion rotation, a being the swirl vector. Hence, eq. ()38|) gives 

Bu; = Su (52) 

The evolution of the anisotropy with time is defined by 

S = -(s2+j(aTo)j(w ))T + S (53) 
j(uJ) = - (Soj(^*o) +j(^ )S )r +i(u ) (54) 

% The operator j stands for the vector product, u x uj = j(u)(uj). 
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which couples the symmetric and the antisymmetric parts of the distortion matrix. The 
swirl magnitude reads 




according to eq. since (3% = 0. Its orientation cannot be determined from the data 
because the above equations describe two distinct kinematics corresponding to ±u7 Q that 
cannot be disentangle. According to eq. (|50j) . if (and only if) the rotation u = then 
the distortion vanishes S = since B is either a symmetric or antisymmetric. In other 
words, a planar distortion has necessarily to account for a rotation. 



4-1. Constant distortion 

Among above solutions which show planar kinematics, let us investigate the (simplest) 
one defined by = 0. In such a case, k Q oc uj and Sou = 0. According to eq. (|53|54|55|) . 
linear calculus shows that the distortion is constant 

S = S , uj = uj (56) 

Such a distortion in the Hubble flows produces a rotating planar velocities field with 
magnitude oc H a~ 2 . In the present case, the model parameters can be easily evaluated 
from data. The observed cosmic velocity fields are partially determined by their radial 
component 

v r = (v, -) — cz, r = rim) = ru = aq, r = ct (57) 
r 

where m, z, u, t stand respectively for the apparent magnitude, the redshift, the line of 
sight, the photon emission date of the galaxy and c the speed of the light. According to 
eq. ()2I9I50|) . the radial velocity of a galaxy located at position f is given by 

H + H^r, Hu = ^-u-Su (58) 

Because trS = 0, it is clear that the sum of three radial velocities v r corresponding to 
galaxies located in the sky toward orthogonal directions and at same distance r provides 
us with the quantity H. Hence, simple algebra shows that the sample average of radial 
velocities within a sphere a radius r is equal to 

(v r ) = H(r) (59) 

Therefore, for motions described by eq. (JTJ, the statistics given in eq. (|59~j) provides us 
with a genuine interpretation of Hubble parameter H = H(t). Hence, according to 
eq. (j9l2QI21|) . one obtains the (generalized) expansion factor 

a(t) = exp f H{t)dt (60) 

Hence, the cosmological parameters can be estimated by fitting the data to the function 

<M*)a , Ko ,c = JP{a) = a 2 H/H , \ Q + k + fi = 1 (61) 
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The component of matrix S can be estimated by substituting H in eq. (|58|). and u is 
obtained from eq. ([56)1 . 

It is clear that the above model is derived in the Galilean reference frame, where 
the Euler-Poisson equations system can be applied. Hence, a non vanishing velocity of 
the observer with respect to this frame an produces a bipolar harmonic signal in the 
distribution of data in the sky, which can be (identified and then) subtracted. 

4-2. Discussion 

The dynamics of a homogenous medium and anisotropic moving under Newton gravity 
was already studied by describing the evolution of an ellipsoid 21J. The current approach 
enables us to identify characteristics of the dynamics of the deformation from isotropic 
Hubble law in a more systematic way by mean of the distorsion matrix. 

At first glance, if the planar anisotropic of the space distribution of galaxies within 
the Local Super Cluster (LSC) is stable then the above solution can be used for 
understanding its cosmic velocity fields, k Q being orthogonal to LSC plane. It is well 
known however that the distribution of galaxies is not so homogenous as that, whereas 
this model describes motions of an homogenous distribution of gravitational sources. 
However, such an approximation level is similar to the one which provides us with the 
observed Hubble law, that is included in this model (S = 0). 

5. Conclusion 

The present anisotropic solution of Euler Poisson equations system generalizes the 
Hubble law and provides us with a better understanding of cosmic velocity fields 
within large scale structures as long as Newton approximation is valid. As a result, 
this generalization of Hubble motion implies necessarily an homogenous distribution of 
gravitational sources, as similarly to Hubble law. Because the chronology identifies to FL 
chronology for a vanishing distortion, this model interprets as a Newton approximation 
of anisotropic cosmological solutions. The motions are characterized by means of a 
constant of motion a. Among them, particular solutions can be easily derived for 
a = 0. They describe all planar distortions, in addition of two classes of planar-axial 
distortions with or without rotation. Among these solutions, the one which ensures a 
planar kinematics is of particular interest because it describes constant (eternal) and 
rotational distortions. This solution can be fully determined from observational data 
except for the orientation of the rotation. The sensible result is that the velocity field 
is not potential. It is interesting to note that this model accounts for motions which 
might be interpreted as due to tidal forces whereas the density is homogeneous. It is an 
alternative to models which assume the presence of gravitational structures similar to 
Great Attractor as origin of a bulk flow. 
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